
rm(list = ls())
## llamando a las librerias que (creemos que) vamos a usar
options("install.lock"=FALSE)
if (!require ("pacman")) install.packages ("pacman")

pacman::p_load (dplyr
                , haven
                , ggplot2
                , ggpubr
                , ggthemes  ## si quieren mas themes
                , readstata13
                , readxl
                , sf
                , tidyverse
                , tidyr
                , units
                , viridis ## paleta de colores Viridis
                , wesanderson## p/usar paleta de colores de Wes Anderson
                , stringr
                , RColorBrewer
                , patchwork
                , Rmisc
                , lfe
                , stargazer
                , AER
                , haven
                , skimr
                , modelsummary
                , terra
                , fixest
                , vtable
                , did
                , cowplot
                , grid)

setwd("/Users/mateovasquez/Dropbox/")

#####

base <- read.dta13("master_violence_landtitle_caracteristicas.dta")

#Treatment 	

base <- base %>% 
  dplyr::mutate(post = ifelse(Año>1995,1,0),
                any_titleXpost = any_title*post,
                any_title2 = any_title,
                post2 = post)


#the mean level of violence in neighboring municipalities (1988–1995)
Neighbors <- base %>% 
  dplyr::filter(Año<=1995) %>% 
  dplyr::group_by(codmpio) %>% 
  dplyr::summarise(spatlagguerrsum100_100_1995 = mean(spatlagguerrsum100_100, na.rm = T))

summary(Neighbors$spatlagguerrsum100_100_1995)

base <- left_join(base,Neighbors,
                  by = "codmpio")

#Controls X post
base <- base %>%
  dplyr::mutate(
    altura_post = altura*post, #ELEVATION
    rainfall_post = rainfall*post, #RAINFALL
    coca_0_1_post = coca_0_1*post, # COCA
    presence_mines1560_post = presence_mines1560*post, #GOLD MINES
    num_encomiendas1560_post = num_encomiendas1560*post, #ENCOMIENDAS
    royalroaddist_post = royalroaddist*post, #distance to royal roads
    distancia_mercado_post = distancia_mercado*post, #nearest market towns (km)
    pop95_post = pop95*post, #municipal population in 1995 
    pct_minority85_post = pct_minority85*post, #share of minority population in 1985 (%)
    conflictos_1901_1917_post = conflictos_1901_1917*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    conflictos_1918_1931_post = conflictos_1918_1931*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    Violencia_48_a_53_post = Violencia_48_a_53*post, #presence of armed combat during La Violencia (1948–1953)
    anucraids1971_78_post = anucraids1971_78*post, #frequency of land invasions between 1971 and 1978
    spatlagguerrsum100_100_1995_post = spatlagguerrsum100_100_1995*post, #the mean level of violence in neighboring municipalities (1988–1995)
    ltotal_plots_rfmd_1960_85_post = ltotal_plots_rfmd_1960_85*post, # total number of plots reformed between 1960 and 1985 to proxy for rural grievances.
    totviolencia85_post = totviolencia85*post,
    any_title_placebo_post = any_title_placebo*post) #ojo: La Violencia (1948 - 1958)

# Applying asinh transformation to each specified variable directly in mutate

base <- base %>%
  mutate(
    as_sumattacks_public = asinh(sumattacks_public),
    as_sumattacks_paramilitary = asinh(sumattacks_paramilitary),
    as_sumattacks_guerrilla = asinh(sumattacks_guerrilla))

# Running the regression with felm, including fixed effects and clustering

base <- base %>%
  mutate(Depto_year = paste(coddepto, Año, sep = "_"))

base_placebo <- base %>% 
  dplyr::filter(any_title_placebo==1 |
                  any_title2==1) %>% 
  dplyr::rename("any_titleXpostXplacebo"="any_titleXpost")

table(base_placebo$any_title2,base_placebo$any_title_placebo)

base_reserve <- base %>% 
  dplyr::filter(any_title_reserve==1 |
                  any_title2==1)%>% 
  dplyr::rename("any_titleXpostXreserve"="any_titleXpost")

table(base_reserve$any_title2,base_reserve$any_title_reserve)

m1 <- felm(as_sumattacks_paramilitary ~
             any_titleXpostXplacebo +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m2 <- felm(as_sumattacks_paramilitary ~
             any_titleXpostXreserve +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)
m3 <- felm(as_sumattacks_public ~
             any_titleXpostXplacebo +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m4 <- felm(as_sumattacks_public ~
             any_titleXpostXreserve +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)
m5 <- felm(as_sumattacks_guerrilla ~
             any_titleXpostXplacebo +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m6 <- felm(as_sumattacks_guerrilla ~
             any_titleXpostXreserve +
             altura_post + 
             rainfall_post + 
             coca_0_1_post + 
             presence_mines1560_post + 
             num_encomiendas1560_post  + 
             royalroaddist_post + 
             distancia_mercado_post + 
             pop95_post + 
             pct_minority85_post + 
             conflictos_1901_1917_post +
             conflictos_1918_1931_post +
             anucraids1971_78_post + 
             totviolencia85_post +
             spatlagguerrsum100_100_1995_post + 
             ltotal_plots_rfmd_1960_85_post  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)

mean_paramilitary_base_reserve <- round(mean(base_reserve$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_reserve <- round(mean(base_reserve$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_reserve <- round(mean(base_reserve$sumattacks_public, na.rm = TRUE), 3)

mean_paramilitary_base_placebo <- round(mean(base_placebo$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_placebo <- round(mean(base_placebo$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_placebo <- round(mean(base_placebo$sumattacks_public, na.rm = TRUE), 3)


stargazer(
  m2,m1,
  m4,m3,
  m6,m5,
  type = "text",
  keep = "any_titleXpost",
  title = "Law 70 Collective Land Titling and Political Violence",
  align = TRUE,
  omit.stat = c("rsq", "ser"),
  covariate.labels = c("Law 70/Resguardo (0/1) × Post (>1995)",
                       "Law 70/Campesino (0/1) × Post (>1995)"),
  dep.var.labels = c("Paramilitary Attacks", 
                     "Police and Army Attacks", 
                     "Guerrilla Attacks"),
  add.lines = list(
    c("Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Municipality F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Controls × Post", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Department × Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Outcome Mean", 
      mean_paramilitary_base_placebo, mean_paramilitary_base_reserve,
      mean_guerrilla_base_placebo, mean_guerrilla_base_reserve, 
      mean_public_base_placebo, mean_public_base_reserve)))

m1w <- felm(as_sumattacks_paramilitary ~
             any_titleXpostXplacebo |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m2w <- felm(as_sumattacks_paramilitary ~
             any_titleXpostXreserve |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)
m3w <- felm(as_sumattacks_public ~
             any_titleXpostXplacebo |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m4w <- felm(as_sumattacks_public ~
             any_titleXpostXreserve |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)
m5w <- felm(as_sumattacks_guerrilla ~
             any_titleXpostXplacebo  |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo)
m6w <- felm(as_sumattacks_guerrilla ~
             any_titleXpostXreserve  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve)

mean_paramilitary_base_reserve <- round(mean(base_reserve$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_reserve <- round(mean(base_reserve$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_reserve <- round(mean(base_reserve$sumattacks_public, na.rm = TRUE), 3)

mean_paramilitary_base_placebo <- round(mean(base_placebo$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_placebo <- round(mean(base_placebo$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_placebo <- round(mean(base_placebo$sumattacks_public, na.rm = TRUE), 3)


stargazer(m1w, m2w, m3w, m4w, m5w, m6w,
  type = "text",
  keep = "any_titleXpost",
  title = "Law 70 Collective Land Titling and Political Violence",
  align = TRUE,
  omit.stat = c("rsq", "ser"),
  dep.var.labels = c("Paramilitary Attacks", 
                     "Police and Army Attacks", 
                     "Guerrilla Attacks"),
  covariate.labels = c("Law 70/Resguardo (0/1) × Post (>1995)",
                       "Law 70/Campesino (0/1) × Post (>1995)"),
  add.lines = list(
    c("Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Municipality F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Controls × Post", "No", "No", "No", "No", "No", "No"),
    c("Department × Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Outcome Mean", 
      mean_paramilitary_base_placebo, mean_paramilitary_base_reserve,
      mean_guerrilla_base_placebo, mean_guerrilla_base_reserve, 
      mean_public_base_placebo, mean_public_base_reserve)))

# Running the regression with felm, including fixed effects and clustering

base_placebo_pacifico <- base %>% 
  dplyr::filter((any_title_placebo==1 |
                  any_title2==1) ) %>% 
  dplyr::rename("any_titleXpostXplacebo"="any_titleXpost")

table(base_placebo$any_title2,base_placebo$any_title_placebo)

base_reserve_pacifico  <- base %>% 
  dplyr::filter((any_title_reserve==1 |
                  any_title2==1))%>% 
  dplyr::rename("any_titleXpostXreserve"="any_titleXpost")

table(base_reserve$any_title2,base_reserve$any_title_reserve)

m1g <- felm(as_sumattacks_paramilitary ~ any_titleXpostXplacebo |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo_pacifico)
m2g <- felm(as_sumattacks_paramilitary ~ any_titleXpostXreserve |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve_pacifico)
m3g <- felm(as_sumattacks_public ~ any_titleXpostXplacebo|Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo_pacifico)
m4g <- felm(as_sumattacks_public ~ any_titleXpostXreserve  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve_pacifico)
m5g <- felm(as_sumattacks_guerrilla ~ any_titleXpostXplacebo  |Depto_year + codmpio + Año | 0 | codmpio, data = base_placebo_pacifico)
m6g <- felm(as_sumattacks_guerrilla ~any_titleXpostXreserve  |Depto_year + codmpio + Año | 0 | codmpio, data = base_reserve_pacifico)

mean_paramilitary_base_reserve <- round(mean(base_reserve_pacifico$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_reserve <- round(mean(base_reserve_pacifico$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_reserve <- round(mean(base_reserve_pacifico$sumattacks_public, na.rm = TRUE), 3)

mean_paramilitary_base_placebo <- round(mean(base_placebo_pacifico$sumattacks_paramilitary, na.rm = TRUE), 3)
mean_guerrilla_base_placebo <- round(mean(base_placebo_pacifico$sumattacks_guerrilla, na.rm = TRUE), 3)
mean_public_base_placebo <- round(mean(base_placebo_pacifico$sumattacks_public, na.rm = TRUE), 3)


stargazer(
  m2g,m1g,
  m4g,m3g,
  m6g,m5g,
  type = "text",
  keep = "any_titleXpost",
  title = "Law 70 Collective Land Titling and Political Violence (Just Pacific Region)",
  align = TRUE,
  omit.stat = c("rsq", "ser"),
  covariate.labels = c("Law 70/Resguardo (0/1) × Post (>1995)",
                       "Law 70/Campesino (0/1) × Post (>1995)"),
  dep.var.labels = c("Paramilitary Attacks", 
                     "Police and Army Attacks", 
                     "Guerrilla Attacks"),
  add.lines = list(
    c("Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Municipality F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Controls × Post", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Department × Year F.E.", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    c("Outcome Mean", 
      mean_paramilitary_base_placebo, mean_paramilitary_base_reserve,
      mean_guerrilla_base_placebo, mean_guerrilla_base_reserve, 
      mean_public_base_placebo, mean_public_base_reserve)))
